home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / multicompare.pro < prev    next >
Text File  |  1997-07-08  |  11KB  |  381 lines
  1. ; $Id: multicompare.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6.  
  7. function testmeans, A
  8. SX = size(A)
  9. if SX(0) eq 1 THEN CNum = 0 else CNum = SX(2) - 1
  10. for i = long(0), CNum do BEGIN
  11. here = where(A(*,i) NE 0,count)
  12. if count NE 2 THEN return,0
  13. ENDFOR
  14. return,1
  15. END
  16.   
  17.  
  18.  
  19.  
  20.  
  21.  
  22.  pro multicompare, X, Contrast, a, ST,One_Way = One, Unequal_One_Way = Unequal,$
  23.             Two_Way = Two, Interactions_Two_Way  = $
  24.             Interactions,LSD = LS, Bonferroni=Bon,Scheffe = Sch, Tukey = Tuk,$
  25.             Missing = M, Block=B, TInteraction = IN
  26. ;+
  27. ; NAME: 
  28. ;    MULTICOMPARE
  29. ;
  30. ; PURPOSE:
  31. ;    Multicompare gives the user access to a variety of procedures for 
  32. ;    making many inferences from a single experimental data set. These
  33. ;    procedures are designed to guard against experimentwise errors 
  34. ;    resulting from the increased probabilty of at least one inferential
  35. ;    error when many inferences are made. 
  36. ;
  37. ; CATEGORY: 
  38. ;    Statistics.
  39. ;
  40. ; CALLING SEQUENCE:
  41. ;    MULTICOMPARE, X, Contrast, A, St
  42. ;
  43. ; INPUTS: 
  44. ;    X:    2 or 3-dimensional array of experimental data values.
  45. ;
  46. ;     Contrast:    An array, dimensioned (CL,C), where:
  47. ;        CL = the number of treatments
  48. ;           = the number of columns of X
  49. ;        and 
  50. ;        C = the number of contrasts or inferences to be tested.
  51. ;        Contrast, B(i,j) = the coefficient of the mean of the ith 
  52. ;        treatment in the jth contrast.
  53. ;
  54. ;    A:    Experimentwise significance level desired. 
  55. ; 
  56. ; OUTPUT:
  57. ;    St:    An array, dimensioned (2,C), where C is the number
  58. ;        of contrasts to test. ST(0,j) = the absolute value
  59. ;        of the test statistic. ST(1,j) = the lower
  60. ;        limit of the rejection region - i.e.,
  61. ;        if ST(0,j) > ST(1,j) reject the null hypothesis for
  62. ;        the jth contrast at the a*100% significance level.
  63. ;
  64. ; KEYWORD PARAMETERS:
  65. ;    Type of design structure. Options are:
  66. ;      ONE_WAY:    If set,  1-way anova.
  67. ;
  68. ; UNEQUAL_ONE_WAY:  If set,  1-way ANOVA with unequal sample sizes.
  69. ;
  70. ;      TWO_WAY:    If set,  2-way ANOVA.
  71. ;
  72. ; INTERACTIONS_TWO_WAY:    If set, 2-way ANOVA with interactions.
  73. ;    One and only one of the above options may be set.
  74. ;
  75. ;    Options for multicomparison testing:
  76. ;
  77. ;    LSD:    Fisher's LSD procedure. LSD is a post-hoc test of any
  78. ;        contrasts of the treatment means.  It should only be used if
  79. ;        the F-test for equal means rejects the null hypothesis at
  80. ;        the A*100% significance level. LSD has an experimentwise
  81. ;        error rate approximately equal to 5%.
  82. ;
  83. ;   BONFERRONI:    Bonferroni's method. This method should be selected
  84. ;        instead of LSD if the F-test is not significant.
  85. ;                                      
  86. ;        
  87. ;      SCHEFFE:    Scheffe's procedure.  Use Sheffe's procedure to make any
  88. ;        number of unplanned comparisons - that is, to data snoop.
  89. ;        Experimentwise error rate is A*100%.              
  90. ;
  91. ;    TUKEY:    Tukey`s procedure. Tukey's procedure is often more sensitive
  92. ;        than Sheffe's but in the general case it requires equal sample
  93. ;        sizes. Pairwise testing of means is allowed with unequal 
  94. ;        sample sizes, but if the disparity is too great, Sheffe's 
  95. ;        method is more sensitive.  The experimentwise error rate, A,
  96. ;        must be between 0 and 0.1.
  97. ;
  98. ;    BLOCK:    a flag which ,if set, signals that comparisons
  99. ;        should be done on block means instead of treatment means.
  100. ;        Alternatively, input transpose(X) for the experimental data
  101. ;        array.
  102. ;
  103. ;      MISSING:    Placeholding value for missing data.  If not set, then assume
  104. ;        no missing data.
  105. ;
  106. ; TINTERACTION:    A flag, if set, to signal that contrasts are for interaction
  107. ;        effects. The user should not set both keywords BLOCK and
  108. ;        TINTERACTION.
  109. ;-
  110. On_Error,2
  111.  
  112. ;  Ascertain design structure
  113.  ONE_WAY = 0 & ONE_WAY_UNEQUAL=1 & TWO_WAY = 2 & TWO_WAY_INTERACTION = 3
  114.  
  115.  keyset = [keyword_set(One), keyword_set(Unequal), $
  116.             keyword_set(Two),keyword_set(Interactions)]
  117.  T  = where(keyset,nkey)
  118.  
  119.  if nkey NE 1 THEN BEGIN
  120.     print, $
  121.          'multicompare - must specify one and only type of '
  122.     print,'               of design structure'
  123.     goto,DONE
  124.  ENDIF
  125.  
  126.  T = T(0)      ; T = selected design structure
  127.  
  128.  
  129. ; Ascertain test type
  130.  LSD = 0 & BONFERRONI = 1 & Scheffe = 2 & TUKEY = 3
  131.  
  132.  keyset = [keyword_set(LS), keyword_set(Bon), $
  133.            keyword_set(Sch),keyword_set(TUK)]
  134.  MT  = where(keyset,nkey)
  135.  
  136.  if nkey NE 1 THEN BEGIN
  137.     print,'multicompare - must specify one and only one type'
  138.    print, '               of multiple comparison test'
  139.    goto,DONE
  140.  ENDIF
  141.  
  142.  MT = MT(0)    ; M = selected test type
  143.  
  144.  
  145. If N_Params(0) LT 3 THEN BEGIN
  146.   print,"multicompare- need more parameters."
  147.   goto, done
  148. ENDIF
  149.  
  150.  
  151. if (a GT 1.0 OR a LT 0) THEN BEGIN
  152.     print,        $
  153.      'multicompare- significance level a must lie between 0 and 1'
  154.     goto, DONE
  155.  ENDIF ELSE if MT EQ TUKEY and a GT .1 THEN BEGIN
  156.       print, "multicompare - experiment-wise error rate must be"  
  157.       print, "               no greater than .1 for Tukey's test"
  158.       goto,done
  159.      ENDIF
  160.  
  161.  X1 = float(X)
  162.  SX = size(X1)
  163.  
  164.  if (SX(0) NE 2 AND T ne TWO_WAY_INTERACTIONS) OR   $
  165.     (SX(0) ne 3 and T eq TWO_WAY_INTERACTIONS) THEN BEGIN
  166.     print,      $
  167.      'multicompare- data array has wrong number of dimensions'
  168.     goto,done
  169.  ENDIF
  170.  
  171. if (KEYWORD_SET(B) NE 0) THEN    $ 
  172.                         ;transpose to test block means
  173.   if (T NE TWO_WAY_INTERACTIONS) THEN      $
  174.        X1 = transpose(X1)        $
  175.   ELSE BEGIN
  176.        SX = size(X1)
  177.        X1 = fltarr(SX(3),SX(2),SX(1))
  178.        for i = 0l,SX(1) -1 do $
  179.          for j = 0l,SX(3) -1 do $
  180.             X1(j,*,i) = float(X(i,*,J))
  181.   ENDELSE
  182.  
  183.  
  184.  SX =size(X1)
  185.  SC = size(Contrast)
  186.  C = SX(1)
  187.  CN = SC(1)
  188.   
  189.  
  190.  
  191.  
  192.  
  193. ; If testing interactions, reform array and perform 1-way testing
  194.  
  195. if KEYWORD_SET(IN) THEN BEGIN
  196.  if KEYWORD_SET(B) THEN BEGIN
  197.   print,   $
  198.     "multicompare - keywords Block and Interaction should not both be set."
  199.   goto,done
  200.  ENDIF 
  201.  
  202.  if (T EQ TWO_WAY_INTERACTIONS) THEN $
  203.    if (SX(0) NE 3) THEN BEGIN
  204.     print, 'multicompare- wrong number of dimensions.' 
  205.     goto,DONE
  206.    ENDIF ELSE BEGIN                ;reform array for one-way anova
  207.    X1 = Fltarr(C*SX(3),SX(2))
  208.    for i = 0,SX(3) -1 DO $
  209.     for j = 0,C-1 DO $
  210.      X1(C * i +j,*) = float(X(j,*,i))
  211.    SX = size(X1)
  212.    C  = SX(1)
  213.    T = ONE_WAY
  214.   ENDELSE     $ 
  215.  ELSE BEGIN
  216.    print,    $
  217.  "multicompare- to compare interaction effects, keyword structure"
  218.   print,"                      must be set to TWO_WAY_INTERACTIONS "
  219.    goto, done
  220.  ENDELSE
  221. ENDIF
  222.  
  223.  
  224.  if ( SC(0) EQ 1) THEN CNUM = 1 else CNUM = SC(2)
  225.  
  226.  if (C NE CN) THEN BEGIN
  227.    print,         $
  228.    'multicompare- data and contrast arrays must have same number of columns'
  229.    goto,DONE
  230.  ENDIF
  231.  
  232.  if testcontrast(Contrast) EQ 0 THEN BEGIN
  233.   print,'multicompare- Contrast array contains equations whose coefficients'
  234.   print,'              do not sum to 0'
  235.   return
  236.  ENDIF
  237.  
  238.  D=1
  239.  if T NE TWO_WAY_INTERACTIONS THEN  R=SX(2) ELSE BEGIN
  240.     R=SX(3) 
  241.     D=SX(2)    
  242.  ENDELSE                                  ;Compute number of rows
  243.  
  244.  
  245.  ; replace missing data with 0
  246.  if ( T EQ ONE_WAY_UNEQUAL and N_ELEMENTS(M) NE 0) THEN BEGIN
  247.  
  248.      Pairwise,X1,M,YR,YC, nogood, good
  249.      CN = where(YC NE 0,NYC)
  250.      RN = where(YR NE 0,NYR)
  251.      if NYC LT 2 OR NYR LT 2 THEN BEGIN
  252.        print,"multicompare, too many missing entries"
  253.       return 
  254.     ENDIF
  255.  
  256.      YC = float(YC)
  257.  
  258.         ; compute sqd dev of columns from their means
  259.      Mean =  (X1 # Replicate(1.,R)) /YC  
  260.      Mean = MEAN # replicate(1.0, R)
  261.      sigma = fltarr(C,R)
  262.      sigma(good) = (X1(good) - Mean(good))^2
  263.      sigma = sigma # Replicate(1.0,R)
  264.  
  265.  ENDIF ELSE    BEGIN
  266.  
  267.          YC = float(R)
  268.         if (T EQ TWO_WAY_INTERACTIONS) THEN BEGIN    
  269.  
  270.           YC = YC * D
  271.           X2 = X1
  272.           X1= reform(X1,C,YC)
  273.   
  274.         ; Compute sum squared deviations of columns from their means 
  275.           Mean =  (X1 # replicate(1.0,YC))/YC
  276.           sigma = (X1 -  MEAN # replicate(1.0, YC))^2 
  277.           sigma = sigma # replicate(1.0,YC)
  278.  
  279.         ;Pre-process three dimensional are by converting to a two
  280.         ; dimensional array whose entries are sums of repetitions
  281.  
  282.         X1 = fltarr(C,R)
  283.         for i = float(0),R-1 Do   $
  284.         for j=0,C-1 DO        $
  285.           X1(j,i) = Total(X2(j,*,i)) 
  286.                           
  287.       ENDIF ELSE BEGIN
  288.        Mean = (X1 # replicate(1.0, YC))/YC
  289.        sigma = (X1 -  MEAN # replicate(1.0, YC))^2 
  290.        sigma = sigma # replicate(1.0,YC)
  291.       ENDELSE
  292.   ENDELSE
  293.  
  294.  
  295.  
  296.  
  297. if ( T EQ ONE_WAY_UNEQUAL AND N_ELEMENTS(M) ne 0) THEN         $
  298.                sigma =sqrt( total(sigma)/(Total(YC)-C))           $
  299. ELSE sigma =sqrt( total(sigma)/(C*YC-C)) 
  300.  
  301.  
  302. Means = (X1#replicate(1.0,R))/YC      ;compute column means
  303. ST = fltarr(2,CNUM)
  304.  
  305. if CNUM GE 2 THEN ST(0,*) = ABS(MEANS # Contrast)  $
  306.  ELSE   ST(0) = ABS(Total(Contrast * Means))
  307.  
  308. CASE MT of
  309.         ; compute test statistic for all examples
  310.         ; based on contrasts
  311.  TUKEY: BEGIN
  312.          If T EQ ONE_WAY_UNEQUAL AND   $
  313.           testmeans(Contrast) EQ 0 THEN BEGIN
  314.           print, 'Multicompare- TUKEY test can only handle'
  315.           print, '              missing data when comparing means'
  316.           goto,done
  317.          ENDIF
  318.  
  319.         if CNUM GE 2 THEN BEGIN 
  320.           if ( T NE ONE_WAY_UNEQUAL OR N_ELEMENTS(M) eq 0) THEN  BEGIN
  321.             Test_Stat =Replicate(1.0,C)   $
  322.                        #ABS(Contrast)/sqrt(YC)   
  323.             ST(1,*) =  sigma * Test_Stat
  324.           ENDIF ELSE BEGIN
  325.  
  326.             for i = 0l,CNUM - 1 DO   $
  327.              ST(1,i) = min (YC(where(Contrast(*,i) NE 0)))
  328.             ST(1,*) =2.0* sigma / sqrt(ST(1,*))
  329.           ENDELSE
  330.        ENDIF ELSE BEGIN
  331.            if T EQ ONE_WAY_UNEQUAL THEN $
  332.             YC1 =  min(YC(where(Contrast ne 0))) $
  333.            else YC1 = YC
  334.            ST(1) = sigma*Total(ABS(Contrast)/sqrt(YC1))
  335.          ENDELSE
  336.   END
  337.  
  338. ELSE: BEGIN
  339.     if CNUM GE 2 THEN BEGIN 
  340.       ;compute test statistic for all examples based on contrasts
  341.     if ( T NE ONE_WAY_UNEQUAL or N_ELEMENTS(M) eq 0) THEN       $
  342.      Test_Stat =Replicate(1.0/YC,C)# Contrast^2   $
  343.     ELSE Test_Stat = 1./YC # Contrast^2
  344.     ST(1,*) = sigma * (sqrt(Test_Stat))
  345.     ENDIF ELSE  $
  346.       ST(1) = sigma*(sqrt(Total(Contrast^2/YC)))
  347.     END
  348. ENDCASE
  349.  
  350.  
  351.  if( N_ELEMENTS(YC) EQ 1) then v = YC*C-C else v = total(YC) -C
  352.  
  353.  case MT of
  354.   
  355.       LSD       : tav = student_t (a/2.,v)
  356.       BONFERRONI:BEGIN
  357.                    tav = student_t(a/(CNUM*2),v)
  358.                    END
  359.       SCHEFFE    : Begin
  360.                      tav =F_test(a,C-1,v)  
  361.                      tav=sqrt((C-1)*tav) 
  362.                      END
  363.        TUKEY:    if C ne 2 THEN $
  364.                   tav = .5 * studrange(1.0 - a,v,C)   $
  365.                 else tav = sqrt(2)*student_t(a/2.,v)/2.
  366.               
  367.        ELSE: BEGIN
  368.             print,"multicompare- unknown Test selected"
  369.             goto,DONE
  370.             END
  371.  ENDCASE
  372.  
  373.  
  374.  
  375. ST(1,*) = tav*ST(1,*)   
  376.  
  377.  
  378. DONE:
  379. RETURN
  380. END   
  381.